Defects dynamics following thermal quenches in square spin-ice 
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I. INTRODUCTION 

Frustrated magnets are classical and quantum systems in 
which the interactions in combination with the lattice struc- 
ture impede the spins to order in an optimal configuration at 
zero temperature 1 1 ]. In classical instances, the local minimi- 
sation of the interaction energy on a frustrated unit gives rise 
to a macroscopic degeneracy of the ground state. This occurs 
in spin-ice samples in which the spin interactions mimic the 



frustration of proton positions in water ice. The theoretical in- 
terest in these systems has been boosted in recent years by the 
artificial design of materials. For instance, regular arrays of 
elongated single-domain ferromagnetic nano-islands arranged 
along the sides of a 2D square lattice, were manufactured. 
The beauty of artificial spin-ice (ASI) is that the state of a sin- 
gle degree of freedom can be directly visualised with magnetic 
force microscopy. Moreover, advances in lithography allow 
great flexibility in their design, and the interaction parameters 
can be precisely controlled by tuning the island length, the 
distance between them, and the height between layers, to se- 
lect the phase into which the system should set in (2H4). One 
of the main goals of the research on artificial spin-ices is to 
develop new materials that could improve the performance of 
data storage and data processing devices [00|. 

One can model spin-ice materials by taking into account 
dipolar interactions I71JT01 or by using a simpler vertex 
model EJHUEl. I n perfect spin-ice samples the total spin 
surrounding a lattice point is constrained to vanish according 
to the two-in/two-out rule, and the vertex model is integrable. 
For the model defined on a square lattice two ordered and a 
critical disordered (spin liquid) phase have been found with 
powerful analytic tools ||T3j[T4l. The model with four-in or 
four-out arrows can also be solved analytically and the criti- 
cal character of the disordered phase is lost in this case [15]. 
Vertices that do not satisfy the two-in/two-out rule and carry 
dipolar moment break integrability and no exact tool exists to 
solve models with them. We recently solved the statics of the 
sixteen-vertex (all possible states of four arrows attached to a 
central site) model with an extension of the Bethe-Peierls or 
cavity method 1 12) and we found intriguing relations between 
this technique and the ones of integrable systems. 

In this work, we will be interested in characterising the sin- 
gle spin flip stochastic dynamics of the sixteen vertex model 
on a two dimensional square lattice with energetically unfa- 
vored defects. Indeed, classical natural frustrated magnets are 
subject to thermal fluctuations, and these can be captured by 
vertex models coupled to a heat-bath (TT] [T6l . In a previ- 
ous Letter we showed that their stochastic dynamics display 
metastability in the disordered phase, coarsening of stripes in 
the ferromagnetic phase, and growth of domains, that can be 
made isotropic, in the antiferromagnetic phase [llj. Here, we 
extend this analysis in several directions to be described in de- 



tail in the main text. We work with system sizes L = O(10 2 ) 
that are of the same order as the ones manufactured with 
lithography techniques when preparing artificial spin-ice sam- 
ples. Finite size effects observed in the simulations are then 
relevant to observations in these materials. 

The manuscript is organised as follows: we first present 
the model and the simulation method in Section [TTJ next, in 
Section 



III 



we discuss the possibility to set the system in 
metastable states with constant number of vertices of each 
kind and an excess of defects following different quenches 
into the PM (S ection [IITA| , a-FM (Section [IITB] ) and c-AF 
(Section IIIQ phases. In Section IV we report our results 
on the coarsening dynamics occurring in the model following 
the before mentioned procedure into the a-FM phase (Sec 
tion 




Figure 1. (Colour online.) The original lattice V, with L 2 vertices, is 
shown in grey. Its medial lattice, V, with 2L 2 spins is shown in red. 



IVB) and c-AF phase (Section [TVCj. Finally, in Sect. |V| The Cartesian system of coordinates of the original lattice (u. 



we present our conclusions. 



II. MODEL AND METHODS 

In this Section we present the model and the numerical 
techniques used in this work. 

A. The sixteen- vertex model 

Conventional vertex models are defined on a finite dimen- 
sional lattice, typically a square one. The degrees of freedom 
(Ising spins, q— valued Potts variables, etc.) sit on the edges of 
the lattice. In our case we use an L x L bidimensional square 
lattice V with unit spacing and periodic boundary conditions 
(PBC). The midpoints of the edges of the lattice V constitute 
another square lattice, the medial lattice V. In the following 
we label the sites of V with (i, j), see Fig.[l| In the model we 
focus on, the degrees of freedom are arrows aligned along the 
edges of the original square lattice V. We can think of them as 
Ising spins Sij = ±1. Without loss of generality, we choose a 
convention such that s = +1 corresponds to an arrow point- 
ing in the right or up direction, depending on the orientation of 
the link; conversely, s = — 1 corresponds to arrows pointing 
down or left. When the spins joining at each vertex of V are 
constrained to satisfy the two-in/two-out rule 021 HI), as for 
the first six vertices in Fig. [2| one has the six-vertex model. 
When the next two vertices in this figure are also allowed 
(four-in and four-out arrows) one has the eight-vertex model. 
Otherwise, the model is generalized to allow for all possible 
vertices with four legs and becomes the sixte en-vertex one, 
that we have already studied in Qj] Q21 033. A charge, q, 
defined as the number of out-going minus the number of in- 
coming arrows, can be attributed to each single vertex config- 
uration l20ll2Tll . Accordingly, in the six- vertex model vertices 
have zero charge and all other vertices breaking the ice-rule 
have a net (positive or negative) charge. 

The energy of each vertex configuration is quantified by 
the Hamiltonian H = J2k=i e k n k, where n^ is the number 
of vertices of type k and e& its energy. We assign a (not 
normalized) Boltzmann weight uj^ = e _/3efc to each of the 
k = 1, . . . , 2 4 four-arrow vertex configurations (note that ujk 



shown in the right bottom end and the (7r/4)-rotated Cartesian system 
of coordinates used to compute parallel and perpendicular correlation 
functions is displayed on the lattice. 

can be greater than one if e& is negative). We set uj\ = co>2 = a, 
co>3 = U4 = b, CJ5 = ujq = c for the ice-rule vertices and 
ujj = 00$ = d, ujq = ... = ujiq = e for the 2-fold and 1-fold 
defects, respectively, ensuring invariance under reversal of all 
arrows (see Fig. [2]). 

6 vertex model 



±± ±± ±± ±± 

a = LUi=UJ2 b=LU3=LU4 C=LU5=UJq d=U>7=U)8 

+ -f + *f + ~h + + 



> 8 vertex model 



16 vertex model 



Figure 2. (Colour online.) The sixteen vertex configurations on the 
2D square lattice and their weights. The first six vertices verify the 
ice-rule. The next pair (2-fold defects) completes the eight-vertex 
model and have charge q = ±2. The last eight vertices (1-fold de- 
fects) have charge q = ±1. 

The static phase diagram of this model was obtained in lfT2ll 
by using the cavity analytic method (Bethe-Peierls approxi- 
mation) and numerical simulations on 2D. We will not repeat 
all the results found in this (long) paper but we simply recall 
that for d ^ and e ^ the model has a complex phase 
diagram with many phases. In particular, we singled out a 
conventional disordered paramagnetic (PM) phase, two ferro- 
magnetic (FM) phases (a and b dominated, respectively), and 
two antiferromagnetic (AF) phases (c and d dominated, re- 
spectively). All transitions are of second order if d ^ and 
e / 0. Using numerical simulations, FM order dominated 
by a vertices was found for a > b + c + d + 3e, close to 
the critical line predicted by the analytic calculations. Using 
the same arguments, AF order dominated by c vertices was 
found for c > a + b + d + 3e and the PM phase was found for 
a, c,d,e < |(a+6+c+<i+3e) (for the precise prediction, that 
we give here only in approximate form, see [ fQ). For small 
defect weight, e, d <C a, 6, c the disordered phase is very close 
to the critical spin-liquid phase of the six- vertex model. The 
closeness to criticality will play an important role in quenches 



into the disordered phase. Henceforth we measure the weights 
in units of c: c — > 1 and a, 6, d, e — > a/c, 6/c, d/c, e/c. 



B. Stochastic dynamics 

We mimic the effect of thermal fluctuations in spin-ice sam- 
ples by coupling the model to an environment and allowing 
for local single spin flips determined by the heat-bath rule. 
Local moves that break the spin-ice rule are not forbidden and 
we therefore allow for thermally-activated creation of defects. 
The dynamics do not conserve any of the various order pa- 
rameters, and are ergodic for both fixed and periodic bound- 
ary conditions. We establish a Monte Carlo algorithm and we 
define the unit of time as a Monte Carlo sweep (MCs). In 
systems with frustration, as the one we are dealing with, com- 
puter time is wasted by the large rejection of blindly proposed 
updates. In order to make the computer time dynamics faster 
we use a rejection-free continuous-time Monte Carlo (MC) 
algorithm [22] . The longest time that we reached with this 
method, once translated in terms of usual MC sweeps, is of 
the order of 10 25 MCs, a scale that is unreachable with usual 
Metropolis algorithms. 

Other kinds of dynamic rules have been used in the liter- 
ature, with different purposes. For instance, a rule that pre- 
serves the ice constraint does not create defects. In order 
to sample the whole phase space on a system with PBC in 
this way one needs to introduce loop updates of any size and 
winding number. Such a dynamics have been studied in the 3- 
colouring model on the hexagonal lattice and leads to glassy 
behavior (23l |24). Another possible local dynamics which 
preserve the ice rules would be to update the system by small 
loops made by four spins around a square plaquette. These dy- 
namics are not ergodic for PBC but they are for the six- vertex 
model with domain-wall BC (DWBC) (25H26). For the spin- 
ice problem, these two possible dynamical models seem quite 
artificial and do not allow us to study defects' motion in the 
way that it is observed to occur in the laboratory. We therefore 
attach to the moves described in the previous paragraph. 

In terms of a reaction-diffusion model the relevant pro- 
cesses taking place during the time-evolution are: 

(2g) + (-2g) -+ (q) + (-g), AE a 2k B T\n(d/e), (1) 
(2q) + (-q) -► (q) + (0), AE ex k B T\n(d/a), (2) 



( q ) + ( q )^(2q) + (0), 

(g) + (-g)->(0) + (0), 



AEock B Tln(e 2 /bd), (3) 
AE ock B Tln(e 2 /ab), (4) 



where the energetic change AE associated with each reaction 
is shown. As AE = E F - E^ with E F = ksThiujp 1 
the energy of the final configuration and Ej = ksThiujJ 1 
the energy of the initial configuration, one has AE = 
ksThi(uji/ujF). In our simulations we will typically use 
a, b ^> e > d as this choice is more relevant experimentally. 
In the first case, eq. ([I]), two defects of type 7 and 8 meet 
to produce two singly (and oppositely) charged defects with 
an energetic gain, if e > d, which depends on the ratio e/d. 
The total density of defects remains constant after this reac- 
tion although their type changes. An example of the second 



case, eq. ([2]), is shown in Fig. [3] a defect of type 7 (charge 
q = 2) meets one of type 14 (charge q = —1) to produce 
a defect of type 10 (charge q = 1) and a spin-ice vertex of 
type 2 with no charge. This corresponds to an energetic gain 
AE < which depends on d/a. Note that the number of sin- 
gle charged defects has not been modified during this process 
but the number of doubly charged defects diminished and so 
did the total number of defects. The reaction in the third line 
represents an initial state made of the 13th vertex (on the left) 
and the 15th vertex (on the right) that turn into a state with 
the 3rd vertex (on the left) and the 7th vertex (on the right) 
by reversing the internal link. The energy variation is then 
AE = /c#Tln(e 2 /6d), and this can be positive or negative 
depending on e/b (smaller than one) vs. e/d (larger than one). 
The fourth reaction is realised, for example, by the reversal 
of the spin on the edge linking vertex number 13 (on the left) 
with vertex number 14 (on the right), leading to vertex number 
3 (on the left) and vertex number 2 (on the right). The energy 
change is then AE — k B T \n(e 2 /ab), a negative quantity for 
our choice of parameters. 

In conclusion, with our choice of parameters, the reactions 
in eqs. ([T}, ^ and Q lead to a decrease in energy while the 
reaction in eq. ^ may be energetically favourable or not de- 
pending on the ratio e 2 /bd. 
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Figure 3. (Colour online.) The reaction in eq. J2|. 



C. Quench dynamics 

We will analyse the system's evolution after an infinitely 
rapid quench from a disordered initial condition into the dis- 
ordered (D), a-ferromagnetic (a-FM) and c-antiferromagnetic 
(c-AF) phases. In practice, we choose a completely disordered 
configuration (T -^oo,a = 6 = d = e = l)asan initial 
condition; such a state is constructed by placing arrows at ran- 
dom on each edge of the square lattice V. If we impose PBC, 
the number of positive and negative charges is identical. We 
subsequently evolve the MC code with parameters that belong 
to the three interesting phases. The system remains globally 
neutral during the evolution, since it is updated by single spin 
flips which cannot create any excess of charge. 

After a quench into the disordered phase the system could 
be expected to equilibrate relatively rapidly; still, it was shown 
in |[TT1l that it remains blocked in metastable states with a finite 
density of defects for long times, if the weight of defects is low 
enough. We will investigate this problem in depth here. We 
will demonstrate that metastability also exists after quenches 
into the ordered phases. Eventually, the interactions between 
the spins, mediated by the choice of vertex weights, creates 
ordered domains of FM or AF kind. The quantitative charac- 
terisation of growth in the ordering processes is given by two 



possibly different growing lengths extracted from correlation 
functions along orthogonal directions || and _L that we identify 
in Fig. [1] 

D. Observables 

The relaxation dynamics of clean lattice systems are usually 
studied in terms of time-dependent macroscopic observables 
averaged over different realisations of the dynamics (thermal 
noise, initial conditions) denoted by (...). In particular, we 
compute the following quantities: 

(i) The density of vertices of each type: 

n a (t) = (n 1 (t)+n 2 (t)) ,n b (t) = (n 3 (t) + n 4 (t)) , (5) 

n c (t) = (n b (t)+ne(t)), (6) 

16 

n d (t) = (n 7 (t)+n 8 (t)) ,n e (t) = ($^*W> > ( 7 ) 

k=9 

n def (t) =n d (t) +n e (t) . (8) 

(ii) The two-times self-correlation function defined by: 

c(t,t w ) = ^ E <%;)W%;)(^)> (9) 

with t > t w . The indices (i, j) denote the coordinates of an 
Ising spin in the medial lattice V (i.e. the vertices of the square 
lattice shown in red in Fig.[T]). 

(iii) The space-time correlation functions. The definition of 
the relevant correlation functions between different points in 
the lattice is not straightforward when we introduce some 
anisotropy in the model (for example, by choosing a > b). 
For convenience, we define a set of correlation functions be- 
tween spins using two different orientations: along the Carte- 
sian axes u x and u y and along the 7r/4-rotated axes u\\ and u± 
(see Fig. [I}. The space-time self correlation functions along 
u\\ and u± are defined as 



G^{r = n,i) = -2 J2 (%i)%i+*0 



(10) 



(hj)ev 



G x (r = n,i) = —^ ^ (S iid) S ii+nd) ) , (11) 



(ij)ev 



where n G N. 



(iv) The growing lengths L"(t) and L^(t) along u\\ and u± 
are extracted numerically from the scaling of the space-time 
correlations: 



(12) 



G^(r,t)c±F Gh J-^- 

V l ii,j 



it)) ' 



III. DEFECT DENSITY 

In this Section we study the density of defects, vertices of 
type e and d after quenches into the different phases. With 
this analysis we investigate the possibility of finding long- 
lived metastable states after dynamic quenches from a fully 
disordered initial condition. 



A. Quench into the PM phase 

In the following, we study the evolution of the model after 
a quench from a random initial condition (a = b = d = 
e = 1) into a different PM state, typically close to the SL 
critical phase (i.e. a = b = 1 and d, e <C 1). In the initial 
configurations defects are common: we are interested here in 
the mechanisms leading to their annihilation. 



1. Equally probable defects, d — e. 

The evolution of the system after a sudden quench into the 
PM phase with d = e was already reported in I'll]. Let 
us recall some useful results: for values of d large enough 
(d = e > 10 ~ 4 ) the density of defects n de f quickly saturates 
to its equilibrium value. For smaller d's, the system gets ar- 
rested into a metastable state with finite and constant density 
of defects n de f = n p for long periods of time. This dynam- 
ical plateau lasts longer for smaller ds, a behaviour reminis- 
cent to what was found in 3D dipolar spin-ice [27]. The time 
regime where the density of defects finally leaves the plateau 
and reaches its equilibrium value, is characterised by a scal- 
ing of the dynamic curves with the characteristic time d~ 2 . 
This scaling strongly suggests that the relevant time scale in 
the system is the typical time needed to create a pair of sin- 
gle defects. From an ice-rule state, the energy change asso- 
ciated with the reaction: (0) + (0) — )■ (q) + (— q) is, e.g., 
AE oc —k B T\n(d 2 /ab), with a = b = 1 in the simula- 
tion. Then, by a simple Arrhenius argument, the typical time 
needed in order to overcome this barrier is ex ex.p(/3AE) giv- 
ing the before mentioned time-scale r ex d~ 2 . 

At a first sight, one could think that the emergence of this 
dynamical plateau in the density of defects is due to the pres- 
ence of doubly charged defects, the ones with q — ±2 and 
weights ujj = ujs = d. Indeed, doubly charged defects 7 and 
8 must decay into two single charged defects in order to be 
able to move. However, the inverse reaction (2q) + (0) — » 
(q) + (q) in eq. fe\ is accompanied by an increase in energy, 
AE ex kBTln(b/d) > when d = e, and it is energetically 
unfavourable. Therefore, <i-vertices get naturally stuck in the 
sample and are very hard to eliminate. This mechanism could 
give a justification for the plateau in the case d = e. 

In real spin-ice realisations, both in 2D and 3D, the energy 
associated to doubly charged defects is much larger than the 
one of single charged ones ej$ ^> eg v .. 5 i6. It is then more 
relevant to experiments to study in detail the effect of d < e 
in the time evolution of the model, and to revisit the influence 
of e and d defects on the development of the plateau. 



2. Single charged defects are more favourable than doubly 
charged ones, d < e. 

We investigate now the dynamical consequences of choos- 
ing different weights for the two kinds of defects. We focus 
on the fate of the dynamical plateau when doubly charged 
defects are rapidly suppressed, with d the smaller weight 
in the model. The inspection of the reaction rates for the 
annihilation-creation of defects suggests to study the two fol- 
lowing cases separately: 

• d < e and d > e 2 : Single charged defects e are slightly 
more favourable than d-defects. However, the decay 
of d-defects into two e-defects following the inverse 
reaction eq. ^ must overcome an energy barrier, as 

AE = k B T\n(bd/e 2 ) > 0. 

• d < e and d < e 2 : Doubly charged defects are very 
unfavourable. The decay of d-defects into two e-defects 
is energetically favoured and occurs spontaneously. 
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Figure 4. (Colour online.) Quench into the PM phase. Decay of the 
density of vertices for a = b = 1 and d < e for a square lattice with 
L — 50 averaged over 500 realisations. The weights of the defects 
are indicated on the figure. The panels in the first row, (a) and (b), 
are for parameters such that d > e 2 , with e 2 /d — 10 -4 in (a) and 



e 2 /d 
ford - 
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J2 



in (b). The panels in the second row, (c) and (d), are 



As shown in Fig. p] (a), for d > e 2 and large enough val- 
ues of d, d/e 2 > 10 , the decay of n e (green data points) and 
rid (blue data points) freeze at a metastable density for around 
five decades in time. The density of e-vertices n e is smaller 
than rid m the plateau regime. Instead, still for d > e 2 but 
for smaller values of d, d/e 2 < 10 2 , d- vertices rapidly disap- 
pear and the plateau is only seen on n e , as shown in panel (b). 
Indeed, after a rapid decay, n e gets frozen into a metastable 
value for a long time before it finally reaches its equilibrium 
value. Hence, one can conclude that the presence of d-defects 
in the system is not responsible for the emergence of the dy- 
namical plateau in the total density of defects ridef- 



The data in panels (c) and (d) in Fig. [4] were obtained for 
d = e 2 . The density of n e remains larger than rid during the 
whole evolution in both cases. Similarly to what was observed 
for e = d EH, the system gets blocked into a metastable 
plateau only for small enough values of e < 10 ~ 4 , and the 
existence of this arrested dynamical regime is not due to an 
excess of d- vertices. The evolution of n e and rid for d < 
e 2 shown in Fig. [4] supports this observation. Although rid 
rapidly vanishes, n e exhibits a dynamical arrest. The value of 
the plateau density can, in principle, depend on the weight of 
the vertices in a complicated manner. We did not study this 
last point in detail here. 

In short, the system can be arrested in long-lived metastable 
states with many defects of type e. 
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Figure 5. (Colour online.) Quench into the PM phase, (a) Time 
evolution of the density of defects, rid and n e , for different boundary 
conditions: periodic boundary conditions (PBC) and fixed boundary 
conditions (FBC). The data were obtained with an average over 500 
realisations of the dynamics with a = b = 1 and d = e=10 -7 for 
a system with linear size L = 50. (b) The time-dependence of the 
total density of defects, ridef = rid + n e , for different system sizes 
with PBC at d — e 2 = 10 -14 , a = b = 1. The linear sizes are given 
in the key and vary by 25 from L = 50 to L — 200. The data were 
obtained after averaging over 300 runs. The initial decay of ridef is 
confronted to po/(l + Ht), with Q a fitting parameter (red dotted 
line), see the text for a discussion. 



3. Boundary conditions 

In order to better understand the emergence of the frozen 
regime we repeated the numerical experiment with fixed 
boundary conditions (FBC): the state of each spin on the 
boundary is kept fixed from the one it had in the (random) ini- 
tial configuration during the simulation. One has to be careful 
when choosing the boundary conditions and make sure that 
these do not induce a polarisation of the sample. Indeed, 
polarised boundary conditions such as the DWBC can have 
dynamical consequences such as the drift of defects (loosely 
speaking, magnetic monopoles). These effects have not been 
studied here. 

In the initial high temperature state, defects of any kind 
populate the system. After the quench, one of the mechanisms 
for relaxation is the annihilation of oppositely charged defects. 
In order to do so, defects have to meet in the appropriate man- 
ner, meaning that the reversal of the spin shared by both of 
them restores the ice rule. In the reaction-diffusion language 
this corresponds to the process (q) + (— q) — >• (0) + (0). Two 



defects of opposite charge ±1 can also meet in the 'wrong' 
way and create, by a single spin-flip, a pair of doubly charged 
defects accordingly to: (q) + (— q) — » (2q) + (— 2q). Starting 
from a completely ordered FM configuration, one can create 
a pair of defects by flipping a string of spins. The latter can 
wind around the lattice by PBC. Then, in order to annihilate 
these pair of defects, one must flip back all the spins belong- 
ing to the string. One can imagine that this kind of extended 
structures could be responsible for the slowing down of the 
dynamics. If so, the evolution of the system with FBC, where 
winding strings are absent, should not present a dynamical 
plateau. As shown in Fig.|5]this is not the case: a metastable 
plateau in the evolution of the density of defects appears with 
FBC as well. This is due to the fact that, in the presence of 
more than a single pair of defects, there is always a way to an- 
nihilate all the defects without going through the boundaries 
of the lattice. In this sense, the dynamics are insensitive to the 
nature of the boundary conditions. 



4. Finite size effects 

As already pointed out in |[TTlL the metastable density of 
defects for d = e depends on the linear size of the system. 
One should then ask, for generic parameters, whether the ob- 
served metastable density is just a finite size effect or not. In 
order to give an answer to this question we simulated systems 
of different sizes under the same conditions that we chose to 
be a = b = 1 and d = e 2 = 10 -14 . The results obtained 
are shown in Fig. [5](b). The height of the plateau, n p , and 
the time spent by the system in this regime decreases with the 
size of the system. However, as the plateau height is subject 
to strong fluctuations, we are not able to predict its precise de- 
pendence on the system size. The data do not show saturation 
at a finite value n p nor length of the plateau and this gives a 
strong indication that this effect is due to the finiteness of the 
samples. This is confirmed by the data obtained after a quench 
into the c-AF phase for different system sizes (see Fig. [8]): in 
this phase the data are less noisy and we found that the plateau 
density n p depends on the size of the system as ex L -2 which 
vanishes in the thermodynamic limit. Having said this, we 
wish to stress that the system sizes of artificial spin-ice sam- 
ples are of the same order as the ones used in our numerics 
and, therefore, blocking effects of the kind here shown are ex- 
pected to exist in those samples as well, e.g. L « 150 in the 
samples studied in (3. 



5. Initial decay 

For d = e, see [11], the initial decay of ridef is well fitted 
by a power law decay 



n def(t) = 



Po 



(i + nt)<* 



(13) 



density agrees with the mean-field reaction-diffusion picture, 
7idef(t) ~ Po /(I + ^)> proposed in [27) for 3D dipolar spin- 
ice. The presence of d-defects modifies this behaviour and 
makes the decay slower [11], i.e. a decreases. This suggests 
that the exponent a depends on d/e 2 and crosses over from 
a = 1 for d/e 2 <C 1 (d < e 2 ) to a < 1 beyond this limit. 



B. Quench into the FM phase 

We now turn to the ordering dynamics following a quench 
from a random initial condition into the FM phase dominated 
by a-vertices (i.e. a>6+l + <i + 3eas discussed in detail 
in [12]). As shown in Fig. [6] a dynamic arrest occurs for small 
defects' weights during the relaxation towards the a-FM phase 
as well. We anticipate that the same kind of behaviour also 
appears when the system is quenched into the c-AF phase (see 
Fig.®. 




with a = 0.78. The evolution of n e for d = e 2 shown in 
Fig. [5] (b) is rather well fitted by a diffusive decay po/(l + 
Qt). Interestingly, when d < e 2 the decay of the defects' 
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Figure 6. (Colour online.) Quench into the FM phase. Time depen- 
dent density of defects after a quench from a — b — d — e— lio 
a = 5, b = 1 and d — e 2 for the different values of e. (a) L — 50 
and e = 10" 3 , 10" 4 , ..., 10" 9 (as shown in the key), (b) L = 100 
for e = 10" 2 (in black), 10" 3 , ..., 10" 10 . For e = 10" 2 the system 
saturates to its equilibrium value shown with a dotted black line, (c) 
Short-time behaviour for e = 10~ 3 , L = 50 and L = 100. The 
decay is confronted to po/(l + fit) (grey dashed line) and the fit 
po/(l + Qt) a with a — 0.59 (blue dashed line), (d) Test of scaling 
with te 2 for L = 50. 

In this section we analyse the relaxation towards the FM 
phase by studying the decay of the defects' density for differ- 
ent values of the external parameters. In Fig. [6] we show the 
evolution of the density of defects after a quench to a = 5, 
b = 1, d = e 2 and different values of e for two different sys- 
tem's sizes L = 50 (a) and L = 100 (b). The data shown 
have been averaged over 10 3 independent realisations of the 
dynamics. For small enough e (e < 10 -3 ) the system gets 
frozen into a metastable state with a finite density of defects. 
Similarly to what is observed after the quench into the PM 
phase [11], the time the system spends in this plateau is longer 
for smaller e. The ordering process following a quench into 



the FM phase is characterised by a time scale r ex e~ 2 in the 
regime in which ridef leaves the plateau. As shown in Fig. [6] 
(d) for L = 50 all the curves collapse into a single curve when 
rescaling the time variable by r. Therefore, the typical time 
associated with the creation of a pair of defects is the relevant 
time scale during this time regime. 



the plateau will disappear in the thermodynamic limit. The 
initial decay is algebraic with a non-trivial power t~ a with 
a = 0.74 independently of the system size. 
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Figure 7. (Colour online.) Quench into the FM phase. Time- 
dependent density of defects after a quench from a—b—d—e—\ 
to b = 1, d = 10" 18 , e = 10" 6 and different values of a. (a) L = 50 
for a = 1, 3, 5, 7, 10, 15 (as shown in the key), (b) L = 100 for 
a = 2, 4, 7, 9. The data have been averaged over 300 independent 
runs. 

The evolution of the density of defects following a quench 
into different points of the a-FM phase is shown in Fig. [7] for 
L = 50 (a) and L = 100 (b) samples. During a short time 
regime (t < 10 MCs) the density of defects decays indepen- 
dently of a. For later times, the decay depends on the value of 
a. In particular, the expected power-law decay n(t) ~ t~ a be- 
comes slower for larger values of a. Therefore, the exponent 
a depends on the weights of the vertices and decreases when 
increasing a. The metas table density of defects increases with 
a and depends on the system size. 

The initial decay of ridef can be fitted by the power-law in 
eq. ( fT3) with a ~ 0.59 over the whole time regime before the 
system reaches the plateau density [as shown in Fig. [6](c)]. 
In the a-FM the decay becomes slower than the diffusive law, 
£ -1 , found in the disordered phase. The algebraic decay does 
not depend on the size of the system as suggested by the data 
shown in Fig. [6] 

A general statement can be made at this point: the value of 
the exponent a decreases - and hence the relaxation becomes 
slower - when going deeper into an ordered phase (of FM or 
AFkind). 



C. Quench into the c-AF phase 

We follow now the evolution of the system after a quench 
from a random initial condition into the AF phase dominated 
by c-vertices (i.e. 1 > a + b + d + 3e |12|). Figure [8] dis- 
plays the temporal dependence of the total density of defects, 
ridef, after such a quench for systems with different linear 
sizes given in the key. The inset shows the linear size depen- 
dence of the plateau height extracted from the data in the main 
part of the figure in a double logarithmic scale. The data points 
are accurately fitted by a L~ 2 dependence that suggests that 
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Figure 8. (Colour online. online.)) Quench into the c-AF phase. 
Temporal decay of the density of defects in the system following 
a quench to the parameters d — e 2 — 10 -14 , a = b = 0.1 and 
c = 1. The data plotted were obtained by using different lattice sizes 
L — 50, ..., 250 and averaging over 300 independent realisations of 
the dynamics. The initial decay is confronted to the power-law decay 
t~ a with a = 0.74. 



D. Conclusion 

In this section, we have shown that the system gets arrested 
into a long-lived frozen state for all kind of quenches as long 
as the weights e and d are small enough (d ^ e < 10 -4 ). 
After inspection of the persistence of the plateau for a large 
range of parameters and different boundary conditions, we 
concluded that the emergence of such dynamical arrest is not 
due to the presence of doubly charged defects and it is not an 
artifact of the periodic boundary conditions. The metastable 
density depends on the system size, in such a way that it might 
disappear in the thermodynamic limit. This assumption is sup- 
ported by the L~ 2 dependence of the metastable density of 
defects after a quench into the c-AF phase (see Fig. [5]). The 
scaling ridef (t) ~ f(te 2 ) shown in Fig.pKd) gives a simple 
interpretation of the long time dynamics: the creation of a pair 
of defects is needed to 'unblock' the evolution and allows the 
system to reach its equilibrium state. Another result reported 
in this section concerns the initial decay of defects right after 
the quenches. We have shown that the density of defects fol- 
lows a power-law decay characterized by an exponent a which 
depends on the vertex weights. In the PM phase, one recov- 
ers the mean field t _1 decay for d ^ e 2 , as reported in 3D 
dipolar spin-ice l27l . For d = e a smaller exponent a = 0.78 
was found (H ], indicating that the relaxation in qualitatively 
slower. After a quench into the FM and AF phases the decay 
of defects becomes slower when choosing parameters deeper 
into the ordered regions of the phase diagram. For a = 5, 
b = 1 and d = e 2 = 10 -6 we found a = 0.59 and we further 
showed that, at fixed value of 6, d, e, this exponent decreases 
when increasing a. 



IV. COARSENING DYNAMICS AND AGEING 

In this Section we analyse the ordering process following 
a quench into the FM and AF phases from a totally random 
initial condition. We choose to work with defect weights sat- 
isfying d = e 2 , differently from what we presented in [11], 
where d = e. The interest is to investigate the role payed by 
each kind of defect in the ordering dynamics, for parameter 
that are closer to the experimental ones in ASI [411281. 



A. Slow relaxation towards the PM phase 

In Fig. [9] we show the decay of the two-time correlation 
function C(t,t w )asa. function of the time difference t—t w for 
different values of t w shown in the key, and working parame- 
ters such that a = b = 1 and d = e 2 (with different choices 
of the defect weights in the two panels, d — e 2 = 10 -4 in 
(a) and d = e 2 = 10 -12 in (b)). Recalling the results shown 
in Fig.[4](c) and (d), metastability (a plateau in the number of 
defects) is not expected in (a) as the value ofd = e 2 is suffi- 
ciently large, while it is expected after times, say, of the order 
of 10 4 MCs in (b) since d = e 2 it is pretty small (e < 10 -4 ). 

One can distinguish different dynamical regimes from these 
curves. For short times, as long as neighbouring defects an- 
nihilate in a few MCs the correlations are time translational 
invariant and close to one. At later times time-translational 
invariance is lost and the system exhibits non- stationary re- 
laxation. The longer the waiting-time t w is, the slower the 
decay, as in an ageing situation. 



In panel (a), when d 



is relatively large, a station- 



ary regime is attained for waiting-times of the order of t w ~ 
3 x 10 5 MCs as demonstrated by the fact that the two curves 
corresponding to the longest t w , beyond this time-scale, fall 
on top of each other and do not depend on t w . As shown in 
the figure, the equilibrium curve follows a stretched exponen- 
tial decay: C eq (t,t w ) = F(t - t w ) = exp[-(t - t w )^/r e ] 
with 7 = 0.81 and r e « 25 x 10 3 MCs. One could expect 
a similar picture away from the 'spin ice' (a = b = 1) curve 
d = e 2 in the disordered phase, with some dependence of the 
parameters t eq , the relaxation time, 7, the stretching exponent 
and r e , the characteristic time. We have not studied these de- 
pendencies in detail. 

The behaviour is different in panel (b), where the weight of 
defects is much smaller. The curves do not reach a stationary 
regime for the waiting times used. Moreover, a long plateau 
is seen in the curve for t w « 10 7 MCs, reminiscent of the 
plateau in the decay of n e shown in Fig. |4](d). The system 
does not evolve during a period of time in between « 10 7 and 
10 9 MCs. In this case an extremely long time-decay is neces- 
sary to reach a complete decorrelation. A steady state regime 
is not reached in the simulation, the correlations continue to 
evolve for all waiting times shown. 



As argued in Sec. Ill the fire sinite of the lattice has very 
strong effects in these systems. For larger system sizes the 
cross-over to metastability will be pushed to longer times, and 
to infinity in the thermodynamic limit. The ever-lasting non- 
stationary relaxation is due to the proximity to the spin-liquid 



critical phase of the six vertex model when the defect weights 
are small enough. In any critical relaxation the system in ques- 
tion will grow equilibrium patches with a time-dependent crit- 
ical length. This length will need an infinite time to reach the 
size of the system if this diverged. 
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Figure 9. (Colour online.) Quench into the FM 
initial configuration. Two-time self-correlation 
a system with L = 50 and data averaged over 
system evolves with a = b = 1 and d — e 2 — 
e 2 = 1CT 12 in (b). Note that in both cases d = 
correlation function in (a) is confronted to exp 
7 = 0.81 andr e = 23014 (solid black line), a 
data very accurately. 
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phase from a random 
function C(t,t w ) for 
500 realizations. The 
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B. Anisotropic FM domain growth 



Starting from a random initial configuration, we quench the 
system by setting a = 5, b = 1 and d = e 2 = 10 -10 at t = 0. 
The equilibrium state corresponding to this set of parameters 
is deep into the FM phase. This choice strongly favours a- 
vertices. The system will then evolve towards its ordered FM 
state by growing domains made of type- 1 and type-2 vertices. 

A first hint into the dynamics of the system is given by the 
time evolution of the density of vertices, n K (t), for each kind 
of vertex, k = a, 6, c, d, e. This is shown in Fig. 10 Dif- 
ferently from what we presented in ifTTTl we distinguish here 
the density of each kind of defect as their weights are not the 
same. The data are accompanied by four configurations that 
illustrate the evolution of the system. One clearly observes the 
growth of anisotropic ordered FM domains. The directions u\\ 
and u± defined in Fig. [T] are parallel and orthogonal to the 
longer domain walls, respectively. 



From inspection of the data plotted in Fig. [10] we can iden- 
tify four different dynamical regimes: 

(I) A short time regime (t < 10 _1 MCs) during which all den- 
sities vary very little (see, e.g., the data in IfTTTl where we gave 
more data-points on this very short regime, for a different set 
of parameters). 

(II) An intermediate regime (t < 10 MCs) characterised by 
the annihilation of a large number of defects which are trans- 
formed into ice-rule vertices by a few single spin-flips. Then 
rid and n e decay (independently of the value of a, as shown 
in Fig. [7]) while n a , rib and n c increase. Quite surprisingly n^ 
and n c increase in the same way in this regime. The typical 
configuration shows no apparent order although the tendency 
to align in a diagonal direction is already visible. 

(III) A slow relaxation regime in which the dominant dynami- 



cal mechanism is the one of growing anisotropic FM domains 
made by type-1 (black) or type-2 (white) vertices. The third 
snapshot illustrates this situation: there is the same number of 
type-1 and type-2 vertices. 

(IV) A much slower regime sets in once a FM domain perco- 
lates in the u\\ direction. This regime is characterised by the 
emergence of very stable FM stripes (t > 10 3 MCs). How- 
ever, the sample is still very far from equilibrium as it needs to 
grow order in the u± direction as well in order to fully equi- 
librate. The percolation of a domain in this latter direction is 
achieved by a extremely slow mechanism that we describe be- 
low. 

(V) The system equilibrates at much longer times (t eq > 10 11 
MCs for this system size). 
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Figure 10. (Colour online.) FM ordering. Upper panel: time evolu- 
tion of the density of vertices with weight a,b,c, d, and e for a = 5, 
b — 1, d — e 2 — 10 -10 , in a system with linear size L — 100. The 
data are averaged over 300 samples. The snapshots are typical con- 
figurations at the instants indicated by the arrows. Black and white 
points are vertices 1 and 2 as defined by the colour code shown in 
Fig.0 

The behaviour of the space-time correlation functions con- 
firm this growth. As shown in Fig. [TT] the correlations along 
the direction parallel to the stripes (b) grow faster than in the 
orthogonal direction (a). The function G" decreases mono- 
tonically and does not vanish at any distance for times larger 
than « 10 3 MCs (regime IV for which FM stripes perco- 
late). Instead, the correlations along u± decrease rapidly at 
small distances and show a minimum at L/(2y/2), then G x 
increases. This is due to the tendency of the system to de- 
velop a modulated structure in time. Because of the periodic 
boundary conditions a stripe is constrained to wind around the 
lattice resulting in a modulated configuration (as illustrated 



by the right-most snapshot in Fig. 10). The growth is highly 
anisotropic since a > b. For b > a correlations along u±_ 
develop faster than along u\\ , forming stripes perpendicular to 



the ones shown in Fig. 10 The relevant parameter characteris- 
ing the anisotropy of the ordering process is the ratio a/b. As 



shown in the insets of Fig.[TT] in the regime where anisotropic 
domains grow (regime III with times t < 10 3 MCs) the cor- 
relation function along both directions depends on space and 



time through the ratio r/t 1 / 2 : 



U/2 



(14) 



which confirms the growing length L±\\(t) ~ t 1 / 2 found 
for d = e ifTTll with a less refined analysis. The master 
curves in the insets of Fig. [ll] are confronted to F- 1 '" (x) = 
exp[— (x/v± i \\) w± >w] where x = r/y/t. The best fits shown 
were obtained with v± = 0.89, w±_ = 1.15 (a) and v\\ = 0.28, 
w\\ = 1.15 (b). Notably, the stretching exponents w are the 
same (within numerical accuracy) but the scales v± t \\ are dif- 
ferent. 
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Figure 11. (Colour online.) Quench into the FM phase. Space-time 
correlation function along the orthogonal (a) and longitudinal (b) di- 
rections for L = 100, a = 5, b = 1, d = e 2 = 10" 10 and different 
times t given in the key. The data have been averaged over 300 runs. 
The insets show G -1 '" as a function of the rescaled variable r/t 1 ^ 2 
for different times in the coarsening regime III. The master curves 
in the insets are confronted to F ± '^(x) = exp[—(x/v± i \\) w± >W)] 
where x = r/y/i. The scaling functions F ± '" shown in thick black 
lines were obtained with v± = 0.89, w± = 1.15 (a) and v\\ = 0.28, 
^H= 1.15(b). 

The dynamical processes leading the dynamics during the 
different dynamical regimes can be understood from the anal- 
ysis of the snapshots: 



(i) As illustrated in Fig. [12j the topology of the model - 
i.e. the fact that a vertex of type-1 cannot have a neighbour of 
type-2 - leads to straight domain walls along the u\\ direction 
made by c- vertices. Instead, the system would develop inter- 
faces between FM states along the u± when quenched into 
the 6-favoured FM phase. Plaquettes of ice-rule vertices, as 
shown in the central panel in Fig. [12| can frequently appear 
along the domain walls. The latter are obtained from the el- 
ementary excitations of the system. These 'loop' fluctuations 
are obtained by sequentially flipping the spins around a pla- 
quette, an operation which preserves the ice-rule (see the first 
and second panel in Fig.[T2|). 

(ii) Domains of the same type are connected by quasi-one- 
dimensional paths made of b- and c- vertices (loop fluctuation 
can eventually be attached to them as well) running through 
a region with the opposite order. The interplay between the 
tendency to order and the local constraint gives rise to these 
structures. They are similar to the ones found in the kinet- 
ically constrained spiral model [29|. In order to further in- 
crease the density of a- vertices and develop the FM order, the 
domain walls and strings of b- and c- vertices have to be elimi- 
nated. The latter disappear first via the following mechanism. 
Curved domains must have 'corners' made of b vertices or de- 
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Figure 12. (Colour online.) Interfaces between FM domains. Lo- 
cal spins on the bonds and vertices are shown following the Colour 
rule established in Fig [2] Left panel: diagonal wall (blue solid line) 
separating two regions with opposite magnetisation both in the x- 
and ^/-directions. Central panel: a 'loop' fluctuation (red solid line) 
on the plaquette highlighted in the left panel. Right panel: a b corner 
vertex cannot be a neighbour of an a-vertex, illustrating the necessity 
of an outgoing string at the corner of an ordered domain. 



C. Isotropic c-AF domain growth 

We now turn to quenches into the AF phase. To start with 
we present results for parameters that favour isotropy, that is 
to say a = b, and defect weights such that e > d. This is the 
case realised in as-grown ASI samples |4, 28]. Next we show 
some snapshots for non-equal FM weights, a ^ b, though still 
in the AF phase to demonstrate that anisotropic AF growth is 
also possible. 



1. c-AF sixteen-vertex model with d ■ 




Figure 13. (Colour online.) Schematic representation of FM stripe 
motion. Vertices on each site are specified. Note that the interface 
is made of antiferromagnetic vertices. Diagonal (red) lines delimit 
domains of opposite magnetisation. Black arrows indicate the spins 
that flip to get the new configuration (represented in blue after the 
flip). 

fects, but b vertices on the corners cannot be surrounded by 
more than two type 1 or 2 vertices (only defects can, giving 
rise to the before mentioned quasi one-dimensional structures, 



as illustrated in the third panel in Fig. 12 ). The string then pro- 
gressively disappears eaten by the attached domains that grow 
from the corner or, alternatively, it is first cut by the creation of 
two defects and the two strands subsequently shrink. Once the 
string has been eliminated one is left with two defects sitting 
on the walls of the now detached domains, that move along 
the interface and eventually annihilate with their anti-partner. 



(iii) Once the domains have grown enough to percolate in 

takes over. 



the u\\ direction the mechanism sketched in Fig. 13 
The formation of stable parallel stripes freezes me dynamics. 
The only way to evolve from a configuration with parallel FM 
stripes is to create a pair of defects by a single spin-flip. After 
the creation of a pair of defects on the domain wall, the diffu- 
sion of the latter along the interface shrink one among the two 
oppositely ordered stripes. This is illustrated by the sequence 



of steps in Fig. [13] The number of steps scales as the length of 
the stripe, hence the time associated to this process diverges 
with the size of the system. 
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Figure 14. (Colour online.) c-AF ordering. Time evolution of the 
density of vertices in a system with L = 100 after a quench to a = 
0.1, b = 0.1, d — e 2 — 10 -10 averaged over 500 runs. Snapshots 
taken at the instants indicated by arrows in the main panel are shown 
(using the colour code defined in Fig.[2j. 



In Fig. 14 we show the evolution of the vertex popula- 



tions, n K , following a quench into the c-AF phase by setting 
a = b = 10 _1 and d = e 2 = 10 -10 . This data are illustrated 
by configurations taken at different instants indicated by ar- 
rows in the main panel. This choice of parameters favours 
c-vertices and is invariant under 7r/ 2 -rotations. Therefore, the 
evolution proceeds by growing isotopic domains of opposite 
staggered magnetisation denoted by m x _l y = ±1 (the stag- 
gered orientation of the spins lying on horizontal and vertical 
edges respectively, see Fig. [I}. 

With a similar analysis to the one used for the FM quench, 
different dynamical regimes can be identified: 

(I) A short transient where all the densities remain roughly 
constant. 

(II) An intermediate regime (t < 10 MCs) with a rapid anni- 
hilation of defects into ice-rule vertices. 

(III) A third regime during which n a = n& decreases in order 
to grow isotropic c-ordered domains. The identity of a and b 
vertex weights implies that n a and n^ are equal. This is ex- 
plicitly shown by the data plotted in Fig. [15] The space-time 
self correlation functions along the || and _L directions are al- 
most identical and the associated growing lengths are, within 
numerical accuracy, t x l 2 . The scaling of the space-time cor- 
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relation function shown in the inset to panel (b) confirms this 
claim. Moreover, the scaling function is again a stretched ex- 
ponential with an stretching exponent w that is very close to 
the one found in the FM quenches. The scale v depends on 
the working parameters. For these set of parameters and sys- 
tem size, regime III is relatively short, it lasts until t ~ 5 10 2 
MCs. 

(IV) Although we have not performed a full statistical anal- 
ysis of this fact, it seems to us that the system freezes when 
an ordered structure winds around the finite size sample for 
periodic boundary conditions (or goes from one border to an- 
other for open boundary conditions). This regime is the one 
in which the plateau in the density of defects, discussed in the 
previous Section, emerges. In the plot in which we show the 
densities of each kind of vertices separately one sees that all 
densities are constant in this regime. The process whereby the 
system leaves this regime will be discussed below. The en- 
trance into this regime will take longer times for larger sizes 
as already discussed in Sec. [in] 

(V) At longer times the system finally reaches equilibrium. 




Figure 15. (Colour online.) Quench into the c-AF phase. Space- 
time correlations after a quench from a random initial condition into 



0.1 



0.1,. 



10 and L — 50 averaged over 



500 runs, (a) Correlations along the _L direction as a function of the 
distance r between sites for different times (shown in the key), (b) 
Correlations along the || direction as a function of r for the same 
times as in (a). The inset shows the scaling of G" as a function of 
r J \ft confronted to g(x) — exp[—(x/v) w )] withx = r j \ft (shown 
in thick black lines). The best fit shown in the inset was obtained 
with v — 0.39 and w — 1.17, that is very close to the value found 
for the stretching exponent in the FM phase. 

A better understanding of the processes involved in the or- 
dering dynamics is reached from the analysis of the snapshots. 

(i) Domain walls are made of a- and b- vertices. Contrarily 
to the FM case, domains of any shape can be constructed with- 
out defects. As shown in the left and central panels in Fig. [16 
horizontal and vertical walls are made by alternating a- and 
b- vertices. Diagonal walls are exclusively made by a- or b- 
vertices depending on their orientation. Therefore, domain 
walls without defects (energetically favoured) form loops of 
spins pointing along the same direction. 

Dynamic domain walls in, say, Ising magnetic models with 
non-conserved order parameter (NCOP) dynamics are curved, 
and at finite times they display a variety of shapes, are rel- 
atively smooth at short length-scales scales and have frac- 
tal properties at long length-scales (3(3 [5T]]. Instead, in the 
sixteen- vertex model with d = e 2 <C 1 c-AF domains have 
tendency to form straight walls made by FM vertices as shown 



in the right most snapshot in Fig. [14] A statistical and geomet- 
ric analysis of the morphology of domains and interfaces and 
their dependence on the parameters of the model remains an 
interesting project, especially if one wishes to confront the 
predictions of this model to images of artificial spin ice sam- 
ples. 

(ii) In the Solid-on-Solid (SOS) representation of the six- 
vertex model 1 32], each domain can be interpreted as a con- 
tour line delimiting regions with different height. The order- 
ing then proceeds by growing or shrinking regions of constant 
height. Note that the Kosterlitz-Thouless phase transition of 
the F-model (a = b) can be mapped onto the roughening tran- 
sition [32]. The ordering dynamics in this phase thus corre- 
spond to flatten the initial rough surface by growing regions 
made by c- vertices. 

(iii) Once isotropic domains are created, one has to elimi- 
nate small domains in order to further increase the density of 
c- vertices of the losing kind and develop the conquering c-AF 



order. Figure 17 illustrates the mechanism taking place. After 
the creation of a pair of defects in a typical time ~ 1/e 2 , their 
motion along the wall shrinks the domain. This is done with- 
out any energy cost and the sequence of steps needed to make 
a domain disappear should scale with its size. The same kind 
of mechanism takes place for horizontal domain walls. 
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Figure 16. (Colour online.) Domain walls between c-AF domains of 
opposite staggered order. The configuration of arrows is shown. We 
use the colour rule defined in Fig. [2] Left panel: diagonal walls in the 
|| direction are made of a-vertices. Central panel: diagonal walls in 
the _L direction are made of 6-vertices. Right panel: horizontal (and 
vertical, by symmetry) walls are made of an alternating chain of a- 
and b- vertices. 
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Figure 17. (Colour online.) Schematic representation of the annihi- 
lation of c-AF domains. Vertices on the walls are represented by the 
colour rule defined in Fig. [2] Blue arrows indicate the spins that have 
been flipped to get the new configuration. As defects diffuse along 
the domain wall the ordered region on the right-bottom side of the 
figure shrinks. 



12 



2. Artificial Spin Ice 

As shown in (T9), one can choose the vertex weights c > 
a = b > e > d in order to model as-grown ASI sam- 
ples (UGH). In these experiments, ferromagnetic islands are 
gradually grown by deposition. During the early stages of 
the preparation the system feels thermal fluctuations and thus 
tries to accommodate to its c-AF ground state. Once the is- 
lands reach a critical size the system blocks into a frozen 
configuration. The question as to whether such frozen con- 
figurations sample the equilibrium distribution has been re- 
cently addressed experimentally [4, 28] and theoretically [19]. 
Using this latter choice of parameters inspired by ASI sam- 
ples: a = exp(— /3e a ), b = exp(— /3e&), c = exp(— /3e c ), 
d = exp(— /3ed) and e = exp(— /3e e ), with e a = e^ = 2/1, 
e c = 2(1 - 2y/2)/l 9 e d = 2(2a/2 + 1)// and e e = (I being 
the lattice constant); the AF domain walls become smoother 
than the ones obtained for d = e 2 = 10 -10 . It is due to the 
presence of more defects lying on the interfaces (see Fig. 6 
in lfT9l ). This kind of domain wall pattern has also been ob- 
served in as-grown artificial spin-ice samples [4]. This is why 
we believe that some of these samples, the ones that are close 
to the phase transition, are out of equilibrium although the 
density of defects is notably well described by the model in 
equilibrium [12] and in simulation studies of a spin-ice model 
system lfT6l . 

Therefore, the shape of the interfaces separating opposite 
AF regions is extremely sensitive to the choice of parameters. 
Setting a slightly above (below) b would produce a preference 
to grow domains in the u\\ (u±) direction. This is illustrated 
by the two snapshots shown in Fig. 18 As shown in |[T2lL the 
equilibrium fluctuations of the model reflect this feature as 
well. The evolution of the system after a quench into the c-AF 
phase would eventually get frozen into an extremely slow re- 
laxation regime, due to the presence of percolating domains 
in the u\\ (u±) direction. The anisotropic dynamical scal- 
ing is difficult to study with our numerical methods in this 
case, and it is not clear whether the same scaling (L\\ ~ t 1 ! 2 ) 
holds along both directions. Experimentally, the anisotropy 
ratio a/b can be tuned by placing islands of different length 
along the horizontal and vertical edges. This work should in- 
cite experimentalists to study this rich far from equilibrium 
behaviour in ASI. 



V. CONCLUSION 

In this paper we presented a thorough study of the relax- 
ational stochastic dynamics of the sixteen vertex model with 
low defect weights after quenches from infinite temperature to 
the three equilibrium phases: disordered, ferromagnetic and 
antiferromagnetic ones. 

In (TT][T2l[l9] we proposed that this celebrated model cap- 
tures many important aspects of two dimensional artificial 
spin ice systems. The interest in these materials lies on the 
fact that they are a possible new technology for storage data 
devices. Our approach to these systems is, clearly, theoretical 
as their fabrication and preparation pose a number of rather 





Figure 18. (Colour online.) Typical configurations of an L = 100 
lattice with a = 0.1, b = 0.01, d = e 2 = 10" 10 (left) and a = 0.1, 
b = 0.001, d — e 2 = 10 -10 (right). The two snapshots have been 
taken at different times during the ordering regime IV in the c-AF 
phase. 

fundamental questions that we can try to give an answer to by 
better understanding the behaviour of the model. 

In the present publication we showed that, for sufficiently 
low defect weight and finite though large system size, the 
single spin-flip stochastic dynamics of the model that mim- 
ics thermal agitation in artificial spin ice, after quenches from 
infinite temperature, present several, very distinct, collective 
dynamic regimes. In short, these regimes are the rapid an- 
nihilation of defects at very short time scales, coarsening of 
nearly critical or phase ordering kind in the paramagnetic or 
ordered phases, respectively, metastability over a very long 
period of time, and the final approach to equilibrium. The 
coarsening process in the paramagnetic phase is controlled by 
the proximity of the critical phase in the absence of defects, 
the so-called spin-liquid phase of the six vertex model. The 
phase ordering kinetics in the ferromagnetic phases is char- 
acterised by anisotropic growth of the ferromagnetic equilib- 
rium states related by spin-reversal. The ordering process in 
the anti-ferromagnetic phase is isotropic for parameters be- 
longing to the so-called F-model (a = b) but it is not isotropic 
for other choices of weights. 

In quenches to the FM phase we found that the end of the 
first coarsening regime is attained when two diagonal stripes 
of opposite ferromagnetic order percolate across the sample. 
Note that these states are equivalent from a thermodynamic 
point of view and no interaction in the model's Hamiltonian 
distinguishes between them. It is the dynamics that devel- 
ops this transient, although very long-lived, anisotropy. The 
crossover from coarsening to metastability clearly depends 
on the size of the sample and the growth law t 1 ! 2 suggests 
a crossover- time t co ~ L 2 when an ordered band goes from 
one system border to another. The fact that the dynamics get 
blocked with two stripes is based on the visual analysis of 
the snapshots and the fact that the perpendicular space-time 
correlation reaches a minimum at a distance of the order of 
L/(2y/2). A more refined analysis along the lines in l33l[34l 
is necessary to clarify the role played by the initial condition 
in the number of stripes in the blocked state. As the time 
needed to let the stripes percolate across the sample diverges 
with the system size, it is clear that metastability and the fur- 
ther approach to equilibrium (regimes IV and V) are pushed 
to infinity in the thermodynamic limit. 
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In quenches to the c-AF phase the entrance in the 
metastable regime is also achieved when the two domains with 
the opposite staggered magnetisation percolate and, therefore, 
the last two regimes are moved to divergent time- scales in the 
thermodynamic limit in this phase as well. AF ordering is 
characterised by a t 1 ! 2 growing length. 

In quenches to the PM phase, the proximity to the critical 
spin-liquid of the six vertex model leads to a relaxation with 
a very long relaxation time that diverges in the limit of zero 
defect weight. Equilibrium patches are hard to visualise in a 
model with many local variables as the single vertices we have 
in this model but, most probably, the cross-over to metastabil- 
ity is also controlled by percolation of one (or more) of these 
regions in this phase. A critical growing length that should 
saturate to a finite value should exist in this regime as well. 

We found stretch exponential scaling functions for the 
space-time correlation functions in the two anisotropic direc- 
tions during FM coarsening and in the isotropic AF coarsen- 
ing. The fits are consistent with the same (within numerical 
accuracy) stretching exponent w but different scale v, i.e. for 
a scaling function F{x) ~ exp[— (x/v) w ). 

Long-lived metastable states exist in the three phases af- 
ter a cross-over time that should diverge with the system size. 
Still, as artificial spin-ice samples are of finite size, we stud- 
ied these states in detail here. The configurations found in 
these time scales, that we called regime IV in the Sections 
concerning the ordered phases, are characterised by densities 
of each kind of vertex that are still far from their equilibrium 
values for the parameters that we investigated in this paper 
(note, however, that agreement between dynamic and static 
densities of vertices can be achieved for other parameters, as 
some of the ones used in [12] with the aim of confronting to 
experimental data). In particular, we found here a finite num- 
ber of defects in the dynamic configurations. Our finite size 
scaling analysis suggests that their density vanish in the infi- 
nite size limit though it is hard to reach a definite conclusion 
in this respect. 



Regime IV is succeeded by a totally different epoch dur- 
ing which complete ordering is eventually reached. Several 
mechanisms for the motion of domain walls in the FM and 
AF phases that constitute the relevant relaxation processes in 
this time regime have been identified and discussed in the text. 

Let us conclude by confronting our approach to other dy- 
namic studies of the same physical problem, reported in the 
literature in the last two or three years. Metastability in the 
defect density of frustrated magnets was discussed in |271 by 
using a reaction-diffusion model in which the defects are rep- 
resented by charges and interactions of Coulomb type are con- 
sidered. In this modelling the 'background configuration' is 
not taken into account. In our model, instead, the structure 
in which the defects are placed and displace is very important 
and determines the crossover time to metastability, the frozen 
nature of the system, the escape time from this frozen state 
and the subsequent final relaxation to equilibrium. In lfT6l , a 
mean-field model of spin-ice, with only short range interac- 
tions, is proposed in order to study the domain wall dynamics 
of ASI in the presence of disorder and magnetic fields, two 
'external' perturbations that we have not taken into account in 
our work. A very different approach has been taken in |[T0l to 
address spin-ice dynamics in 2D. In this work the authors in- 
troduce a model of continuous magnetic moments with dipo- 
lar interactions and study its stochastic dynamics. The strong 
anisotropy imposed along the edges of the underlying square 
lattice makes the connection with ASI samples possible. From 
a more general perspective on 2D frustrated magnetism, nu- 
merical studies of the constrained dynamics of a simple clas- 
sical lattice model have shown the existence of a dynamical 
arrest analogous to the one discussed in the text (24). In this 
work the defects breaking a local constraint (analogous to the 
ice-rules) are strictly forbidden contrarily to what we do in our 
paper. Our work presents a new approach that allows to deal 
with thermal fluctuations breaking a hard constraint, and its 
consequences in the out-of-equilibrium relaxation of geomet- 
rically frustrated magnets. 
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